Lecture_5
Doynikov A
February 17, 2016
Сезонные ARMA модели
При помощи общей модели авторегресии-скользящего среднего \[φ(B)x_{t}=θ(B)\epsilon_{t}\] можно вполне адекватно описывать достаточно широкий класс временных рядов. Но если задаться построением наиболее экономичной (с точки зрения количества оцениваемых параметров) модели, то в некоторых случаях небольшая ее модификация позволяет значительно уменьшить число параметров модели. В первую очередь это касается так называемых сезонных временных рядов, в которых наблюдается определенное сходство в поведении наблюдений, разделенных по времени на фиксированный пероид времени s. Пусть во временном ряде \(z_{t}\) нам удалось связать наблюдения, отстоящие друг от друга на время s в модель \[Φ_{P}(B^{s})∇_{s}^{D}z_{t}=Θ_{Q}(B^{s})a_{t} (*)\]
Где \(Φ_{P}(B^{s})\) и \(Θ_{Q}(B^{s})\) полиномы степени \(P,Q\) соответственно вида \[Φ_{P}(B^{s})=1-Φ_1B^{s}-...-Φ_{P}B^{Ps}\] \[Θ_{Q}(B^{s})=1-Θ_1B^{s}-...-Θ_{Q}B^{Qs}\] которые удовлетворяют условиям стационарности и обратимости, a \[∇_{s}^{D}=(1-B^{s})^{D}\] оператор сезонной разности порядка \(D\)
Ошибки \(a_{t}\) здесь не обязательно должны быть некоррелированными, так как наличие связи между переменными через период времени s не означает, что нет связи между соседними значениями ряда \(z_{t}\) и \(z_{t-1}\) или переменными, отстоящими друг от друга на два, три и т.д. моментами времени. Чтобы учесть и эту связь мы вводим еще одну модель \[φ(B)∇^{d}a_{t}=θ(B)\epsilon_{t}\]
где \(\epsilon_{t}\)- уже белый шум, а \(φ(B)\) и \(θ(B)\) полиномы по \(B\) степени \(p\) и \(q\) соответсвенно, также удовлетворяющие условиям стационарности и обратимости. Подставляя последнее выражение в (*) получим окончательную общую модель \[φ(B)Φ_{P}(B^{s})∇^{d}∇_{s}^{D}z_{t}=Θ_{Q}(B^{s})θ(B)\epsilon_{t} \] (**)
которую и называют сезонной моделью авторегресии проинтегрированного скользящего среднего и обозначают \(SARMA(p,q,d)(P,Q,D)\)
Таким образом мы задаем мультипликативную модель сезонности.
Идентификация сезонных моделей
Сезоннные модели скользяшего среднего
Для идентификация сезонных моделей также широко используется АКФ и ЧАКФ. Ранее мы вводили понятие производящей функции автокорреляций. Напомню это
\[\gamma(z)= \sum_{k=-\infty}^{\infty}c_kz^k\]
Для выяснения характеров их поведения ACF и PACF нам поможет тот факт, что производящая функция автоковариаций процесса заданного мультипликативно равна произведению производящих функций его компонент. Т.е. если компоненты \[φ(B)∇^{d}a_{t}=θ(B)\epsilon_{t}\] \[\Phi_{P}(B^{s})\Delta^{D}z_t=\Theta_{Q}(B^s)a_t\] имеют производящие функции АКФ \(γ(B),Γ(B)\) соответсвенно, то производящая функция автоковариаций процесса \(z_{t}\) равна \(γ(B)Γ(B)\)
Проще всего найти выражение для модели сезонного скользящего среднего. \[x_t=\Theta_Q(B^s)\theta_q(B)\epsilon_t\]. Для такой модели производящие функции автокорреляций равны \[γ(B)=\sigma^2_{\epsilon}\theta_q(B)\theta(B^{-1})\],\[\Gamma(B)=\sigma^2_{a} \Theta_Q(B^s)\Theta_Q(B^{-s})\] Cледовательно для всего процесса она будет \[C(B)=\sigma^2_{\epsilon}\sigma^2_{a}\theta_q(B)\Theta_Q(B^s)\theta_q(B^{-1})\Theta_Q(B^{-s})\]
Пример SARMA(0,1,0)(0,1,0)
В этой модели \(\theta_q(B)=1 - \theta_1B\),\(\Theta_Q(B^s)=1 - \Theta_1B^s\) поэтому \[c(B)= \sigma^2_{\epsilon}\sigma^2_{a}(1 - \theta_1B)(1 - \theta_1B^{-1})(1 - \Theta_1B^s)(1 - \Theta_1B^{-s})\] Перемножая полиномы, получим, что отличны от нуля будут коэффициенты при степени \(0,1,,...,s-1,s,s+1\)
library(polynom)
library(fArma)## Loading required package: timeDate## Loading required package: timeSeries## Loading required package: fBasics## ## Rmetrics Package fBasics## Analysing Markets and calculating Basic Statistics## Copyright (C) 2005-2014 Rmetrics Association Zurich## Educational Software for Financial Engineering and Computational Science## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.## https://www.rmetrics.org --- Mail to: info@rmetrics.org season <- 7
theta1 <- 0.6
Theta1 <- 0.8
ma <- c(1,theta1)
sma <- c(1,rep(0,season-1),Theta1)
ma <- as.polynomial(ma)
sma <- as.polynomial(sma)
ma## 1 + 0.6*x sma## 1 + 0.8*x^7 sma_ma <- ma *sma
sma_ma## 1 + 0.6*x + 0.8*x^7 + 0.48*x^8 ma_vec <- sma_ma[2:length(sma_ma)]
ma_vec## [1] 0.60 0.00 0.00 0.00 0.00 0.00 0.80 0.48Оценим ACF по модели
model <- list(ma=ma_vec)
armaTrueacf(model=model,lag.max = 30,type = "correlation") 
## lag acf
## 0 0 1.0000000
## 1 1 0.4411765
## 2 2 0.0000000
## 3 3 0.0000000
## 4 4 0.0000000
## 5 5 0.0000000
## 6 6 0.2152080
## 7 7 0.4878049
## 8 8 0.2152080
## 9 9 0.0000000
## 10 10 0.0000000
## 11 11 0.0000000
## 12 12 0.0000000
## 13 13 0.0000000
## 14 14 0.0000000
## 15 15 0.0000000
## 16 16 0.0000000
## 17 17 0.0000000
## 18 18 0.0000000
## 19 19 0.0000000
## 20 20 0.0000000
## 21 21 0.0000000
## 22 22 0.0000000
## 23 23 0.0000000
## 24 24 0.0000000
## 25 25 0.0000000
## 26 26 0.0000000
## 27 27 0.0000000
## 28 28 0.0000000
## 29 29 0.0000000
## 30 30 0.0000000Смоделируем сезонный ряд
sma <- armaSim(model = list(ma=ma_vec),n = 300)
matplot(sma,type = "l",col = "blue",main = "Simulated Seasonal Time Series",lwd = 2)
Оценим ACF по ряду
acf(sma,lag.max=30,lwd = 3)
Для сезонного скользящего среднего произвольного порядка SARMA(0,Q,0)(0,q,0) c периодом сезонности \(s\) нетрудно убедиться, что отличны от нуля будут тольк значения ACF на задержках \[0,1,...,q,...,s-q,...,s-1,s,s+1,...,s+q,....,sQ-q,...,sQ-1,sQ,sQ+1,...,sQ+q\]
Для PACF без доказательсва примем на веру тот факт, что PACF процесса сезонного сколящего среднего экспоненциально затухает.
Сезоннные модели авторегрессии
Здесь нам опять помогут уравнения Юла-Уокера из прошлой лекции. Из них немедленно следует, что PACF процесса SARMA(p,0)(P,0) сезонной авторегресссии с периодом сезонности \(s\) вида \[(1-\Phi_1B^s+\Phi_2B^{2S}+...+\Phi_PB^{Ps})x_t=a_t\] \[(1-\phi_1B-\phi_2B^2-...-\phi_pB^P)a_t=\epsilon_t\] будет отлична от нуля на задержках \[0,1,...,p,...,s-p,...,s-1,s,s+1,...,s+p,....,sP-p,...,sP-1,sP,sP+1,...,sP+p\]
Для ACF аналогично доказывается, как и в Лекции 3,только более громоздко, что ACF представляет собой экспоненциально затухающие косинусоиды.
Пример SARMA(1,0,0)(1,0,0)
Модель имеет вид \[x_t=\Phi_1x_{t-s} + a_t\] \[a_t= \phi_1a_{t-1}+\epsilon_t\]
library(polynom)
library(fArma)
season <- 7
phi1 <- 0.4
Phi1 <- -0.3
ar <- c(1,phi1)
sar <- c(1,rep(0,season-1),Phi1)
ar <- as.polynomial(ar)
sar <- as.polynomial(sar)
ar## 1 + 0.4*x sar## 1 - 0.3*x^7 sar_ar <- ar *sar
sar_ar## 1 + 0.4*x - 0.3*x^7 - 0.12*x^8 ar_vec <- sar_ar[2:length(sar_ar)]
ar_vec## [1] 0.40 0.00 0.00 0.00 0.00 0.00 -0.30 -0.12Оценим ACF,PACF по модели
modelSAR <- list(ar=ar_vec)
armaTrueacf(model=modelSAR,lag.max = 30,type = "both") 
## lag acf
## 0 0 1.0000000000
## 1 1 0.5210623842
## 2 2 0.2615446164
## 3 3 0.1121635265
## 4 4 0.0100146006
## 5 5 -0.0879171679
## 6 6 -0.2228709364
## 7 7 -0.4516758607
## 8 8 -0.4569890595
## 9 9 -0.3237864948
## 10 10 -0.1945490099
## 11 11 -0.0942836073
## 12 12 -0.0125400446
## 13 13 0.0723953232
## 14 14 0.1912053998
## 15 15 0.2677799811
## 16 16 0.2590866280
## 17 17 0.2008537335
## 18 18 0.1319724568
## 19 19 0.0678650290
## 20 20 0.0069322200
## 21 21 -0.0632761707
## 22 22 -0.1285891106
## 23 23 -0.1612952304
## 24 24 -0.1558646076
## 25 25 -0.1260400281
## 26 26 -0.0866122147
## 27 27 -0.0448683554
## 28 28 0.0002036427
## 29 29 0.0462513307
## 30 30 0.0823197947Смоделируем сезонный ряд
sarma <- armaSim(model = modelSAR,n = 300)
matplot(sarma,type = "l",col = "blue",main = "Simulated Seasonal Time Series",lwd = 2)
Оценим ACF,PACF по ряду
acf(sarma,lag.max=40,lwd = 3)
pacf(sarma,lag.max=40,lwd = 3)
Помочь освоить идентификацию сезонных моделей поможет следующий фрейм
Оценка сезонных моделей
При оценке сезонных моделей используются те же методы, что и для несезонных моделей рядов. Об этом была Лекция 3. Ведь сезонные ряды это просто экономная (с точки зрения числа параметоров) форма записи обычной ARMA модели. Это методы моментов (Юла-Уокера) для сезонной авторегрессии, метод наименьших квадратов и максимального правдоподобия для произвольных сезонных моделей.
Осталось только продемострировать идентификацию и оценку параметров реального экономического сезонного ряда:Импорт в млрд. долларов Российской Федерации.
Информация о динамике объемов импорта и экспорта отражается в отчете о внешнеторговом балансе государства. Внешнеторговый баланс - разница экспорта и импорта товаров. Если экспорт превышает импорт, то образуется положительное сальдо внешнеторгового баланса или внешнеторговый профицит. Если импорт превышает экспорт, то возникает отрицательное сальдо внешнеторгового баланса или внешнеторговый дефицит. В России экспорт традиционно превышает импорт за счет масштабных экспортных поставок газа, нефти и металлов. Соответственно, сальдо торгового баланса РФ сильно зависимо от мировых цен на данные сырьевые группы товаров.
Читаем данные из файла
import <- read.csv("c:/ShurD/Alection/import.csv", header = TRUE)
head(import)## Data Month Year Fact Forecast X Previous Currency
## 1 19.02.1997 ßíâ 97 4.730 NA ìëðä.$ NA
## 2 19.03.1997 Ôåâ 97 4.965 NA 4.73 ìëðä.$ NA
## 3 23.04.1997 Ìàð 97 5.545 NA 4.965 ìëðä.$ NA
## 4 21.05.1997 Àïð 97 6.062 NA 5.545 ìëðä.$ NA
## 5 18.06.1997 Ìàé 97 5.411 NA 6.062 ìëðä.$ NA
## 6 23.07.1997 Èþí 97 5.694 NA 5.411 ìëðä.$ NAimportFact <- ts(import$Fact ,start= c(1997,1),frequency =12)
plot(importFact ,type = "l", col = "blue",lwd = 2,main = "Import (mlrd dollar).Russia")
Нестационарность явно присутствует в ряде. Возьмем обычную разность \[ y_t= x_t-x_{t-1}\]
y <- diff(importFact)
plot(y ,type = "l", col = "blue",lwd = 2,main = "Import (mlrd dollar).Russia.Differences")
Заметна нестационарность в сезонности с периодом 12. Возьмем сезонную разность \[z_t= y_t-y_{t-12}\]
z <- diff(y,lag = 12)
plot(z ,type = "l", col = "blue",lwd = 2,main = "Import (mlrd dollar).Russia.Seasonal Differences")
Нестационарность по дисперсии (или как сейчас принято говорить по волатильности). Но модели таких рядов мы пока не изучали. Поэтому возьмем выделим подряд с 2006 года, после которого волатильность на вид постоянна.
zz <- z[100:216]
zz <- ts(zz,start= c(2006,5),frequency = 12 )
plot(zz ,type = "l", col = "blue",lwd = 2,main = "Import (mlrd dollar).Russia.Seasonal Differences")
Оценим ACF,PACF
zz <- zz[1:116]
par(mfrow = c(1, 2), cex = 0.9)
acf(zz, lwd = 5,main = "ACF",col = "blue",lag.max=40)
pacf(zz,lwd = 5,main = "PACF",col = "blue",lag.max=40)
Более выразительна PACF. Попробуем оценить модель сезонной авторегрессии. Период сезонности очевидно 12. Порядок сезонной авторегрессии 1. Порядок несезонной авторегрессии 3 или 4, возьмем 4. Оцениваем эту модель методом максимального правдоподобия.
modelimport <- arima(zz,order = c(4,0,0),seasonal = list(order= c(1,0,0),period = 12 ),method = "ML")
modelimport$coef## ar1 ar2 ar3 ar4 sar1 intercept
## -0.21644967 0.03676760 0.27994636 0.13599941 -0.37447931 -0.07129735Насколько значимы оцененные параметры можно косвенно судить по ковариационной матрице оценок
modelimport$var.coef## ar1 ar2 ar3 ar4
## ar1 0.0084915716 0.0014197888 -3.786230e-04 -0.0020564972
## ar2 0.0014197888 0.0082429298 1.385544e-03 -0.0004030656
## ar3 -0.0003786230 0.0013855445 8.158319e-03 0.0015422504
## ar4 -0.0020564972 -0.0004030656 1.542250e-03 0.0084329584
## sar1 0.0006446262 -0.0009488153 3.667948e-05 0.0009792428
## intercept 0.0001462261 0.0003468278 3.671008e-04 0.0003519086
## sar1 intercept
## ar1 6.446262e-04 0.0001462261
## ar2 -9.488153e-04 0.0003468278
## ar3 3.667948e-05 0.0003671008
## ar4 9.792428e-04 0.0003519086
## sar1 7.568049e-03 0.0003589610
## intercept 3.589610e-04 0.0262124644P-Значения для проверки значимости вызывать отдельно или вычислять самостоятельно. Так как оценки ассимптотически нормальные, дисперсии ошибок стоят на диагонале матрицы ковариаций. Поэтому для получения р-значений для проверки значимости делаем следующую команду
(1-pnorm(abs(modelimport$coef)/sqrt(diag(modelimport$var.coef))))## ar1 ar2 ar3 ar4 sar1
## 9.414705e-03 3.427492e-01 9.696386e-04 6.930692e-02 8.363238e-06
## intercept
## 3.298339e-01Другой способ такой проверки это метод coeftest из библиотеки lmtest
library("lmtest")## Warning: package 'lmtest' was built under R version 3.2.4## Loading required package: zoo##
## Attaching package: 'zoo'## The following object is masked from 'package:timeSeries':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericcoeftest(modelimport)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.216450 0.092150 -2.3489 0.018829 *
## ar2 0.036768 0.090791 0.4050 0.685498
## ar3 0.279946 0.090323 3.0994 0.001939 **
## ar4 0.135999 0.091831 1.4810 0.138614
## sar1 -0.374479 0.086995 -4.3046 1.673e-05 ***
## intercept -0.071297 0.161903 -0.4404 0.659668
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Остатки после исключения модели и усредненная сумма квадратов остатков
par(mfrow = c(1, 1), cex = 0.9)
modelimport$sigma2## [1] 3.248453importresid <- modelimport$residuals
plot(importresid ,type = "l", col = "blue",lwd = 2,main = "Residuals")
проведем тест Льюнг-Бокса остатков и постром qq-график
Box.test(importresid, lag = 16, type = "Ljung-Box", fitdf = 2)##
## Box-Ljung test
##
## data: importresid
## X-squared = 20.122, df = 14, p-value = 0.1263 qqnorm(importresid)
qqline(importresid)
Остатки явно не нормально распределены, имеют очень тяжелые хвосты Построим ACF,PACF остатков
par(mfrow = c(1, 2), cex = 0.9)
pacf(importresid,lwd = 5,lag.max = 40, col = "blue")
acf(importresid,lwd = 5,lag.max = 40, col = "blue")
Похоже модель выбрана не лучшим образом. Попробуем другую модель. Порядок сезонного скользящего среднего 1, а не сезонного 3. Повторим оценку уже другой модели.
modelimport2 <- arima(zz,order = c(0,0,3),seasonal = list(order= c(0,0,1),period = 12 ),method = "ML")
coeftest(modelimport2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.258841 0.091417 -2.8314 0.004634 **
## ma2 0.050746 0.087772 0.5782 0.563159
## ma3 0.218558 0.086069 2.5393 0.011107 *
## sma1 -0.665206 0.117588 -5.6571 1.54e-08 ***
## intercept -0.076175 0.067536 -1.1279 0.259359
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow = c(1, 1), cex = 0.9)
modelimport2$sigma2## [1] 2.820544importresid2 <- modelimport2$residuals
plot(importresid2 ,type = "l", col = "blue",lwd = 2,main = "Residuals")
Box.test(importresid2, lag = 16, type = "Ljung-Box", fitdf = 2)##
## Box-Ljung test
##
## data: importresid2
## X-squared = 13.874, df = 14, p-value = 0.4592 qqnorm(importresid2)
qqline(importresid2)
par(mfrow = c(1, 2), cex = 0.9)
pacf(importresid2,lwd = 5,lag.max = 40, col = "blue")
acf(importresid2,lwd = 5,lag.max = 40, col = "blue")
Ну вот и так мы оценили достаточно адекватную модель SARMA(3,0)(0,1) остатки гипотезе, что это белый шум не противоречат, но их нормальность явно под вопросом. Тест Шапиро (см Лекцию 2.) нормальности это подтверждает, точнее гипотеза о нормальности остатков отвергается.
shapiro.test(importresid2)##
## Shapiro-Wilk normality test
##
## data: importresid2
## W = 0.95513, p-value = 0.0006707Сравнение моделей. Информационный критерий Хироцуго Акаике
Для сравнения результатов оценки по нескольким моделям и выбора лучшей из них Х.Акаике (1974) был предложен следующий критерий. В общем случае AIC: \[\mathit{AIC} = 2k - 2\ln(L)\,\] где \(k\) — число параметров в статистической модели, и \(L\) — максимизированное значение функции правдоподобия модели.
Будем полагать, что ошибки модели нормально и независимо распределены. Пусть \(n\) — число наблюдений и RSS- \[\mathit{RSS} = \sum_{i=1}^n \hat{\varepsilon}_i^2\], остаточная сумма квадратов. Далее мы предполагаем, что дисперсия ошибок модели неизвестна, но одинакова для всех них. Следовательно: \[\mathit{AIC}=2k + n[\ln(2\pi \mathit{RSS}/n) + 1]\,\]. В случае сравнения моделей на выборках одинаковой длины, выражение можно упростить, выкидывая члены зависящие только от n: \[\mathit{AIC}=2k + n[\ln(\mathit{RSS})]\,\]. Таким образом, критерий не только вознаграждает за качество приближения, но и штрафует за использование излишнего количества параметров модели. Считается, что наилучшей будет модель с наименьшим значением критерия AIC. Сравним значения критерия для двух моделей, которые мы оценили выше.
modelimport$aic## [1] 481.9704 modelimport2$aic## [1] 468.7346Что потверждает наши предположения, о том, что вторая модель предпочтительней.Стоит отметить, что абсолютное значение AIC не имеет смысла — он указывает только на относительный порядок сравниваемых моделей.